Rotons in interacting ultracold Bose gases 
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In three dimensions, non-interacting bosons undergo Bose-Einstein condensation at a critical 
temperature, Tc, which is slightly shifted by ATc, if the particles interact. We calculate the excitation 
spectrum of interacting Bose-systems, *He and ®^Rb, and show that a roton minimum emerges in the 
spectrum above a threshold value of the gas parameter. We provide a general theoretical argument 
for why the roton minimum and the maximal upward critical temperature shift are related. We also 
suggest two experimental avenues to observe rotons in condensates. These results, based upon a 
Path-Integral Monte-Carlo approach, provide a microscopic explanation of the shift in the critical 
temperature and also show that a roton minimum does emerge in the excitation spectrum of particles 
with a structureless, short-range, two-body interaction. 
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Appearance of a roton minimum [l[ in the excitation 
spectrum is an essential part of the phenomenology of 
superfluidity. Bogoliubov derived Q the excitation spec- 
trum of a weakly interacting boson system for small mo- 
menta 
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where eo{k) denotes the free-particle energy, n = N/V is 
the particle density, g measures the interaction strength, 
and m is the mass of a particle. The two-body interaction 
is approximated by a repulsive (5-potential with strength 
determined by the s-wave scattering length, a. 

The importance of interactions on the macroscopic 
scale is characterised by the gas parameter, 7 = na"^, 
which is proportional to the ratio of the volume occu- 
pied by the particles compared to that available to them. 
Bose-Einstein condensates in ultra-cold gases are usually 
in the weakly interacting limit and away from the 
Feschbach resonance experiments could only explore the 
dilute limit, 7 ~ 10~^ [J]. However the use of Feshbach 
resonances in ®^Rb has facilitated the creation of con- 
densates with tunable values of 7 up to approximately 
8 X 10~^ 0. In condensates of ^Li, values of 7 up to 

50 have been achieved @, although with non- uniform 
density distribution due to trapping. We note here that 
7 is not at all small, ^ 0.2, in the case of liquid "'He 

Interactions also induce a shift in the condensation 
critical temperature, Tc, which was first reported for the 
^He-Vycor system [7| confirming earlier theoretical pre- 
dictions 
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where ri°^ is the critical temperature of an ideal Bose- 
gas with the same density and C is a positive dimension- 
less constant. Different theoretical approaches aimed at 
determining this dimensionless constant resulted in sig- 



nificant discrepancy for some time, with consensus finally 
provided by Monte Carlo simulations [31 ■ 

For small 7, the Bogoliubov spectrum ([T]) accurately 
describes the excitation spectrum, while for larger 7, such 
as in liquid ^He, it deteriorates and a roton minimum is 
observed Path integral Monte Carlo calculations 



111 . Il6j and an experiment 0] have shown that equation 
becomes invalid for 7 > 10~^. For even higher 7. the 
system freezes and the roton minimum goes soft. 

In fact, ATc, has a maximum when 7 « 10^^ as can 
be seen in Fig. [T] Since the onset of condensation is 
associated with the energies of available states, we might 
postulate the peak in ATc to be accompanied by the 
transition to the roton regime as 7 increases. Below we 
substantiate this claim. 

Contrary to earlier studies in which the condensate was 
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FIG. 1. The critical temperature shift ATc, is plotted as a 
function of 7 as derived by Griiter et al. [ill]. Note the differ- 
ent scales on the ordinate. Reproduced through the courtesy 
of the authors of Ref. [ij . The shaded area covers the range 
where the maximum shift is located. 
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rotating fast [17| or particles had laser-induced dipole 
moments [isj . we show that the roton minimum may 
also occur in a homogeneous system where the two-body 
interaction is short-range and structureless. This fact 
may give hope for indirectly observe rotons in Bose- 
Einstein condensates utilising already existing experi- 
mental aparatus. Moreover, our simulations imply that 
the excitation of rotons is also connected to the shift of 
the critical temperature of condensation scrutinized in 
earlier studies. 

Theoretical model: We assume that our system can 
be described using only two-body interaction potentials 
between the constituents and can be modelled by the 
hard-sphere potential with diameter a. We have varied a, 
while keeping n constant, to sample different values of 7. 
The systems were homogeneous and periodic boundary 
condition have been applied. 

Numerical details: We have carried out path integral 
Monte Carlo (PIMC) simulations ^ for both ''He and 
^^Rb in order to determine the energy spectrum, e{k). 
The configurations of the system are sampled via the 
density matrix. Although the exact density matrix is 
only available at high temperature, the low temperature 
density matrix can be constructed by convolution of a 
sufficient number of high temperature density matrices. 
This iteration forms the basis of PIMC simulations. 

In our simulations, we used 128 (''He) or 64 (®''Rb) 
particles with 64 density matrices and ensured that in- 
creasing the number of density matrices did not alter 
the final result. Number densities and temperatures of 



(^''Rb) are utilised. The density for *^Rb is at the lower 
end of those reported in experiments [1, [i^l, while the 
^He density is similar for comparison. Simulations with 
realistic liquid ^He density have also been carried out 
with similar conclusion presented below. 

Analysis: The essential object of our analysis is the 
static pair correlation function, g{r), which is propor- 
tional to the probability of finding two atoms at a 
distance r from each other. PIMC methods are well 
suited for the calculation of g{r) from which the static 
structure factor can also be acquired via ^(A;) = 1 + 
nr I - 1) sin(A:r) dr [U. Finally, S{k) is linked to 
the excitation spectrum by the following relation (22. 23| 



S(k) = -— coth -f^ 
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providing an upper bound on e{k). 

Due to the finite size of the system, g{r) can only be 
calculated over a limited range (repeated use of the pe- 
riodic boundary condition would lead to unphysical spa- 
tial correlation). One may assume that g{r > Tc) = 1 
for some cut-off Vc, however, this approach introduces a 
non-continuous step in g{r) at Tc and leads to unaccept- 
able error in S{k) below ^ r~^ . The truncation error can 
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FIG. 2. The real parts of the first three lowest complex poles 
are shown. The inset exhibits all the zeros over a range of the 
complex plane for different 7. The leftmost curve belongs to 
7 — 0.095, while the rightmost to 7 = 1.91. As 7 increases the 
zeros move towards the imaginary axis and all zeros become 
purely imaginary for 7 = 1.91. 



be suppressed by varying the cut-off, Tc, and taking the 
average of the values of S{k) obtained this way. 

An alternative is to fit a trial function to g{r). The 
trial function used below has been derived from the 
Ornstein-Zernike equation with the Percus-Yevick (PY) 
closure [l^l providing the analytical, but implicit solu- 
tion gpY{r) = C~\Gpy{z)}, where stands for the 
inverse Laplace-transform and Gpyiz) is known explic- 



25l |26[. The function, Gpy{z), has a pole at the 
origin, zq = 0, and infinitely many distinct conjugate 
pairs of poles, Z£ — K£ + iki {I =1, 2, . . . ). The distri- 
bution of these zeros completely determines gpY(f), and 
thus all thermodynamical quantities of the system, for a 
fixed value of 7. Figure[2]depicts the real part of the three 
lowest lying zeros, while the inset shows all the zeros over 
the complex plane in the region [-12, 0] x [— IOOtt, IOOtt]. 
For moderate values of 7 only z\ contributes significantly 
27| and dominates the asymptotic behaviour as 



<7PY(f) = 1 H cos (k\r -f b\)e~ 
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In Fig. [2] the real parts of the first three lowest lying 
poles are shown. One may notice that k\ approaches zero 
rapidly, and therefore the corresponding oscillatory con- 
tribution to g[r) is only weakly damped. The imaginary 
part ki (smaller inset) determine the angular wavenum- 
ber of the oscillatory contribution. The smallest imagi- 
nary part, /ci, falls between tt and Stt providing a charac- 
teristic wavenumber, i.e. a broad qualitative estimate on 
where the roton minimum may occur, i.e. between 7r/a 
and 37r/a. 

In order to check our PIMC results and also to smooth 
the simulation data of g(r) we determine the fitting pa- 
rameters Ci, fci, i5i, and k,\ for each value of 7. Although 
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Ki and ki could be determined from Gpy{z), we treat 
them as free parameters, and compare their fitted values 
to those the pole-structure suggests. The comparison 
provides a test of the accuracy of our PIMC calculation. 
Satisfactory agreement between the fitted function and 
raw data can be seen in Fig. [31 The advantage of using 
the fitted gpvir) is that it can be extrapolated to larger 
system sizes and its derivative may be easily calculated. 
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FIG. 3. The radial distribution function (x) and the fitted 
first order PY approximation (dash-dotted line) are depicted 
for 7 = 0.26. The statistical error is within the size of the 
crosses. Inset shows the corresponding static structure fac- 
tor. The vertical dashed lines indicate the range for the first 
maximum predicted by the PY approximation. 

The key result of our work, the excitation spectra of 
^He and ^'^Rb can be seen in Fig. The maximum in 
ATc occurs at 7 w 0.01, while wo find that the roton mini- 
mum appears at 7 « 0.2. Although the onset of the roton 
minimum does not appear to directly coincide with the 
maximum of ATc, the minimum is preceded by another 
qualitative change, namely the development of a point 
of inflection in e{k). By taking the second derivative of 
e{k), as shown in the inset of Fig. HJ^a), we can determine 
the approximate value of 7 for which the spectrum de- 
velops this inflection point. The data indicate that this 
occurs approximately at 7 « 0.014 which falls exactly 
into the region where the critical temperature reaches its 
maximum. 

The connection between the development of the inflec- 
tion point and the maximum of ATc may be illuminated 
by the following argument. At temperature T the occu- 
pancy of a state with energy £{k) is determined by the 
Bosc-Einstein distribution, f{T,e). Bose-Einstein con- 
densation occurs if the ground state is macroscopically 
occupied, i.e. the number of particles in any excited 
states, Nc^{T), is saturated 
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FIG. 4. The excitation spectra, are shown for ^He (a) 

and *^Rb (b). The values of gas parameter, 7, are 0.262, 
0.176, 0.111, 0.064, 0.008 for *He and 0.373, 0.176, 0.064, 
0.014 for *^Rb. The inset depicts the second derivative of 
e{k) for different values of 7. The ordinate is in units of 
A = Y? I2ra. Data for 7 = 0.014 nearly reaches zero, thereby 
signaling the emergence of an inflection point in 



It is tacitly assumed that the chemical potential has 
reached its maximum value. The critical temperature, 
Tc, is thus determined by the excitation spectrum via 
the density of states, Gie). If G(e) increases for ther- 
mally available states, then the temperature must be de- 
creased to reduce /(T, e) and preserve the validity of in- 
equality (O. The critical temperature should therefore 
decrease as well. The density of states corresponds to 
|c?A:/de|, which is the reciprocal of the slope of the exci- 
tation spectrum. 

For weakly interacting bosons, Bogoliubov's result ([T]) 
indicates that increasing 7 increases the slope of e(fc), 
thus G{e) must decrease around /c « 0. We can therefore 
conclude that the critical temperature of a weakly inter- 
acting Bose-gas ATc must be positive. As 7 increases, 
e(fc) develops an inflection point and becomes convex in 
a certain region. The decrease of the slope in this region 
opens an abundance of excited states for the particles. 



4 



iTidxon peaK 


movnn / 

1 1 1 ClA^J 1 1 / 

maximum / 






roto n peak 


/ roton 


/ minimum 


dos 1 


/ k 


FIG. 5. Qualitative picture of the excitation spectrum (right) 
and the density of states (left). As e(fc) develops a maxi- 



mum (maxon) and a minimum (roton) the density of states 
diverges. 



manifested as a peak in Gie). These excited states can 
be populated at lower temperature, thus reducing the 
population of the ground state. Therefore the system 
has to be further cooled for macroscopical occupation of 
the ground state. Consequently, the critical temperature 
must decrease, resulting in ATc decreasing also. As 7 
increases further, the inflection gives rise to the roton 
minimum in the spectrum for non-vanishing momenta. 
Although the magnitude of the derivative starts to in- 
crease along the sides of the minimum, the density of 
states continues to increase, as there are now three sets 
of k states which contribute to the density of states for a 
given energy (See Fig. [S]). 

Recently, a similar conjecture has been made in the 
context of two-dimensional dipolar systems 28| . The in- 
teraction parameter corresponding to our 7 is the dipole 
coupling D. As D increases the critical temperature de- 
velopes a maximum. The authors mentioned the appar- 
ent coincidence of this maximum and the onset of the ro- 
ton minimum, although they have not investigated this 
in detail. We note that their system and ours are dis- 
similar in two important aspects: dimensionality and 
the nature of interaction. Dimensionality is significant 
in that it influences which phases appear in the system, 
e.g. Berezinksii-Kosterlitz-Thouless phase (iol - fsH rather 
than BEC. However the nature of the interactions -we 
feel- is more fundamental. The two-body interaction in 
our analysis is spherically symmetric, structureless and 
short-range, while their dipole interaction has strong an- 
gular dependence and long range influence. It is far from 
obvious a priori that hard core bosons, therefore, should 
establish a roton structure in the excitation spectrum at 
all, let alone that its appearance should correspond to 
the shift of the critical temperature. 

Experimental verification of the predicted onset of the 
roton minimum may be possible by utilising a Feshbach 
resonance to tune the interaction strength and then us- 
ing Bragg spectroscopy to probe the excitation spectrum. 
The values of 7 achieved by Papp et. al. ^ approach 



the lower end of values for which we predict the inflec- 
tion point in e{k) to occur, while Bragg spectroscopy has 
previously been used to measure the excitation spectrum 
of a weakly interacting BEC [H] . 

As a second scenario, quantum evaporation could be 
adapted for pancake shaped Bose-Einstein condensates. 
This technique proved to be a successful experimental 
method in the case of liquid ^He [s^ . An excitation hav- 
ing enough energy is capable of ejecting an atom from 
the condensate. However, phonons and rotons carry dif- 
ferent momenta at the same energy, therefore the ejected 
atoms have different angular distribution depending on 
which type of excitations they interacted with. Therefore 
detecting the angular distribution of ejected atoms after 
exciting the condensate in a controlled manner, could 
prove the existence of rotons in BECs. 

In this Letter, we have shown that the maximum ob- 
served in ATc of an interacting BEC is related to the 
appearance of an inflection point in e{k). As 7 increases 
the inflection point signals the appearance of a roton min- 
imum, characteristic of the excitation spectrum of e.g. 
superfluid ^He. We have also provided a physical argu- 
ment as to why this happens and how it could be observed 
experimentally. 
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